Il modello SIR permette di studiare la diffusione di una malattia in una rete sociale. La popolazione viene divisa in 3 classi:
S: il numero di individui sani, suscettibili di essere infettati;
I: il numero di individui infetti;
R: il numero di individui rimessi, guariti o morti (in entrambi i casi non ci si può reinfettare).
N definisce la somma di tutti gli individui, quindi la popolazione totale, che è costante.
Si considera inoltre una probabilità beta di contagio e una probabilità gamma di passare da infetto a rimesso.
In modello SIR può essere implementato in un network sociale.
Si suppone che:
1: non vi sia tempo di incubazione;
2: il contagio avvenga, con una probabilità Beta, per contatto diretto;
3: ogni individuo malato abbia una probabilita di guarigione costante per unita di tempo.
Inoltre si presume che ogni individuo abbia la stessa possibilità, per unità di tempo, di entrare in contatto con ogni altro; ogni persona si mescola e si incontra completamente a caso in questo modello.
E' possibile eventualmente implementare diversi tipi di network randomici che simulano parzialmente la rete sociale. Un esempio potrebbe essere il network Erdős-Rényi in cui ogni nodo ha una determinata probabilità di avere un link.
Nel network, inoltre, ogni nodo rappresenta una persona mentre gli edges tra i nodi sono i link tra le persone. Per essere contagiata una persona di tipo S deve avere, necessariamente, un edge con un infetto. In ogni iterazione i vicini di un infetto sono contagiati con una certa probabilità (Beta), mentre gli infetti diventano rimessi con una probabilità gamma. Diversamente dal modello classico con le equazioni, la struttura del network influenza il propagarsi della malattia.
In questo programma ho utilizzato inoltre la matrice di adiacenza e non l'edge list (che sarebbe più efficiente), perchè ho voluto considerare network soprattutto densi e con un numero di persone non maggiore di 100, in quanto, in questo caso, la rete si vedrebbe in modo confuso a causa della risoluzione e della dimensione limitata delle immagini.
import matplotlib.pyplot as plt
import numpy as np
import networkx as nx
from networkx import community
import random
import matplotlib.patches as mpatches
#definisco la funzione che plotta il network
def my_net_draw(G,colormap,s_num,i_num,r_num,numPic):
images.append(nx.draw(G, pos, node_color=colormap, with_labels=True))#questo fa fare le immagini
# Legenda
s_List.append(s_num)
i_List.append(i_num)
r_List.append(r_num)
blue_patch = mpatches.Patch(color='blue', label='SUSCETTIBILI = %.d ' %(s_num))
red_patch = mpatches.Patch(color='red', label='INFETTI = %.d ' %(i_num))
green_patch = mpatches.Patch(color='green', label='RIMESSI = %.d ' %(r_num))
plt.legend(handles=[blue_patch,red_patch,green_patch])
if(len(images)>0):
filename = f'Graph{numPic}.PNG' #f serve a fare crescere numpic, f-string
plt.savefig(filename)
plt.show(block=False)
numPic+=1
filenames.append(filename)
def my_net_draw_finale(G,colormap,s_num,i_num,r_num,numPic):
images.append(nx.draw(G, pos, node_color=colormap, with_labels=True)) # immagini
blue_patch = mpatches.Patch(color='blue', label='SUSCETTIBILI = %.d ' %(s_num))
red_patch = mpatches.Patch(color='red', label='INFETTI = %.d ' %(i_num))
green_patch = mpatches.Patch(color='green', label='RIMESSI = %.d ' %(r_num))
plt.legend(handles=[blue_patch,red_patch,green_patch])
# Assegna biunivocamente un colore in colormap ad ogni nodo (in base allo status: S, I o R )
def my_check_status(status,colormap):
global s_num,i_num,r_num
if status[i] == 0:
colormap.append('blue')
s_num+=1
if status[i] == 1:
colormap.append('red')
i_num+=1
if status[i] == 2:
colormap.append('green')
r_num+=1
return (s_num,i_num,r_num)
#definisco le costanti e le liste principali.
Susceptible = 0
Infectious = 1
Recovered = 2
status = [] # lista che conterrà 0,1 o 2
pop = [] # qui verrà salvato ogni nodo numerato da 0 a n
adj_matr = [] # edges
networkxedges = [] # List of networkx friendly edges
images = [] # imaggini dei network
#richiedo i parametri
scelto = int(input('premere:\n"0" per inserire manualmente i parametri della malattia \n"1" per usare dei paramteri preimpostati \nScelto: '))
if scelto==0:
numnodi = int(input('Inserire il numero di soggetti con cui avviare la simulazione: '))
Beta = float(input('Inserire la probabilità, un float nell"intervallo chiuso 0-1, di contagio: ' ))
gamma = float(input('Inserire la probabilità, un float nell"intervallo chiuso 0-1, di diventare recovered: ' ))
I0 = int(input('Quanti sono gli infetti all inizio?'))
elif scelto==1:
numnodi = 50
Beta, gamma = 0.2, 0.01
I0 = 1
else:
print('errore')
for h in range(0, numnodi):
pop.append(h)
#ER
scelto=int(input('Scegli quale network implementare tra:\n"1" Erdos-Renyi\n"2" Completo \n"3" numero personalizzato di edges \nScelto:'))
if scelto==1:
print("Inserire la probabilità di avere un link tra i nodi")
P = float(input(':'))
G = nx.Graph()
G.add_nodes_from(range(0, numnodi))
adj_matr=np.zeros(numnodi**2).reshape(numnodi,numnodi)#
for i in G.nodes():
for j in G.nodes():
if (i < j):
R = random.random()
if (R < P):
G.add_edge(i, j)
if(G.has_edge(i,j)==True):
adj_matr[i][j]=adj_matr[j][i]=1
networkxedges.append((j,i))
elif scelto==3 :
#completo
print('Hai scelto numero personalizzato di edges ')
possedges = int(numnodi * (numnodi - 1) * 0.5) # per completo max
print('Numero totale massimo di edges possibili', possedges, ".")
numedge = int(input('Inserire il numero di edges: '))
for h in range(len(pop)):
row=[]
for j in range(len(pop)):
row.append(0)
adj_matr.append(row)
for k in range(0, numedge): #num totale
check = 0
while check == 0:
ranrow=random.randint(0, len(pop)-1)
rancolumn=random.randint(0, len(pop)-1)
if adj_matr[ranrow][rancolumn]!=1 and ranrow!=rancolumn: # edges è la matrice di adiacenza con netwedg vedo quali collego
adj_matr[ranrow][rancolumn]=random.randint(0,1)
adj_matr[rancolumn][ranrow]=adj_matr[ranrow][rancolumn] # perchè è simmetrica
networkxedges.append((rancolumn,ranrow)) # vedo cosa collego. Se S <-> I deve esserci edge
check = 1
G = nx.Graph()
G.add_nodes_from(pop)
G.add_edges_from(networkxedges)
elif scelto==2:
print('hai scelto completo')
P = 1
G = nx.Graph()
G.add_nodes_from(range(0, numnodi))
adj_matr=np.zeros(numnodi**2).reshape(numnodi,numnodi)#
for i in G.nodes():
for j in G.nodes():
if (i < j):
R = random.random()
if (R < P):
G.add_edge(i, j)
if(G.has_edge(i,j)==True):
adj_matr[i][j]=adj_matr[j][i]=1
networkxedges.append((j,i))
G = nx.Graph()
G.add_nodes_from(pop)
G.add_edges_from(networkxedges)
pos = nx.spring_layout(G, iterations=200)
else:
print('Errore')
seed=int(input('Inserire seed: '))
pos = nx.spring_layout(G, seed=seed)
#inizio inserendo tutti come scuscettibili
for j in range(len(pop)):
status.append(Susceptible) # mette pop volte 0
#per ogni I0, seleziono un S0 a caso e se non è già passato s I0 lo faccio passare a I0
for n in range(I0):
check2 = 0
while check2 == 0:
starter = random.randint(0, len(pop) - 1)
if status[starter]!=Infectious: #aggiungo a caso gli infetti (n volte 1)
status[starter] = Infectious
check2 = 1
#inizializzo il grafico
pos = nx.spring_layout(G, iterations=200)
numPic = 1
filenames = []
s_List = []
i_List = []
r_List = []
while Infectious in status: #cioè finchè ho 1
s_num=0
i_num=0
r_num=0
infectable = []
colormap = []
conta=[]
for i in range(len(pop)):
my_check_status(status,colormap)
S = 0
I = 0
R = 0
for g in range(len(pop)):
for u in range(0, g+1):
if adj_matr[g][u]==Infectious: # 1 è sia link in adj matr che inf
if (status[g]==Infectious and status[u]==Susceptible) or (status[g]==Susceptible and status[u]==Infectious) :
infectable.append(u)
copia = status.copy()
for j in range(len(status)):
if status[j] == Infectious:
rand = random.random()
if rand < (gamma):
status[j] = Recovered
# INFECTING
for x in range(len(infectable)):
ran = random.random()
if ran < Beta:
status[infectable[x]] = Infectious
if copia != status:
my_net_draw(G,colormap,s_num,i_num,r_num,numPic)
for nnn in range(len(s_List)):
conta.append(nnn)
nnn=nnn+1
# QUESTA PARTE DEVE STARE FUORI DAL CICLO PERCHE' NON SOFFISFA LA CONDIZIONE SECONDO CUI DEVONO ESSERE PRESENTI INFETTI
# ALLA FINE, INFATTI, NON DEVONO ESSERCI INFETTI (o non sarebbe la fine).
colormap = [] # map of colors
s_num=0
i_num=0
r_num=0
for i in range(len(pop)): #solo per lo stato finale
if status[i] == 0:
colormap.append('blue') #s
s_num+=1
if status[i] == 1:
colormap.append('red') #i
i_num+=1
if status[i] == 2:
colormap.append('green')#r
r_num+=1
#images.append(nx.draw(G, pos, node_color=colormap, with_labels=True))
my_net_draw_finale(G,colormap,s_num,i_num,r_num,numPic)
filename = f'Graph{numPic}.png'
plt.savefig(filename)
plt.show(block=False)
filenames.append(filename)
plt.scatter(conta,s_List,color='blue')
plt.scatter(conta,i_List,color='red')
plt.scatter(conta,r_List,color='green')
plt.plot(conta,s_List)
plt.plot(conta,i_List)
plt.plot(conta,r_List)
plt.show()
premere: "0" per inserire manualmente i parametri della malattia "1" per usare dei paramteri preimpostati Scelto: 0 Inserire il numero di soggetti con cui avviare la simulazione: 40 Inserire la probabilità, un float nell"intervallo chiuso 0-1, di contagio: 0.3 Inserire la probabilità, un float nell"intervallo chiuso 0-1, di diventare recovered: 0.1 Quanti sono gli infetti all inizio?1 Scegli quale network implementare tra: "1" Erdos-Renyi "2" Completo "3" numero personalizzato di edges Scelto:1 Inserire la probabilità di avere un link tra i nodi :0.3 Inserire seed: 12345
Il modello sir può essere visualizzato anche in altri modi:
Se si considerano le equazioni che caratterizzano questo modello è possibile automatizzare la risoluzione del sistema di equazioni differenziali e poi stamparle in un grafico time-dependent.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation
Questa parte di codice sfrutta la funzione 'odeint' di scipy per risolvere le equazioni differenziali del SIR e poi plotta il grafico tipico.
scelto = int(input('premere:\n "0" per inserire manualmente i parametri della malattia \n "1" per usare dei paramtri preimpostati \n '))
if scelto==0:
N = int(input('inserire il numero totale della popolazione: '))
I0 = int(input('inserire il numero iniziale di infetti: '))
R0 = int(input('inserire il numero iniziale di recovered: '))
S0 = N - I0 - R0
beta = float(input('inserire tasso di contagio (float tra 0 e 1 ): '))
gamma = float(input('inserire il tasso di guarigione (float tra 0 e 1 ):'))
elif scelto==1:
N = 1000
I0, R0 = 1, 0
S0 = N - I0 - R0
beta, gamma = 0.2, 1./10
else:
print('errore')
t = np.linspace(0, 200, 200) #genero 200 numeri tra 0 e 200
def deriv(y, t, N, beta, gamma):
S, I, R = y
dSdt = -beta * S * I / N
dIdt = beta * S * I / N - gamma * I
dRdt = gamma * I
return dSdt, dIdt, dRdt
y0 = S0, I0, R0 # Condizioni iniziali
ret = odeint(deriv, y0, t, args=(N, beta, gamma)) # dalla documentazine: scipy.integrate.odeint¶: Integrate a system of ordinary differential equations.
S, I, R = ret.T
# Plot
fig = plt.figure(facecolor='w')
ax = fig.add_subplot( facecolor='#dddddd', axisbelow=True)
# funzioni trovate
ax.plot(t, S, 'b', alpha=0.8, lw=2, label='Suscettibili')
ax.plot(t, I, 'r', alpha=0.8, lw=2, label='Infetti')
ax.plot(t, R, 'g', alpha=0.8, lw=2, label='Rimessi')
# Assi del grafico
ax.set_xlabel('Tempo (Day)')
ax.set_ylabel('num pop (su 1000)')
ax.set_ylim(0,N)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(visible=True, which='major', color='white', linewidth=2, linestyle='--')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
ax.spines[spine].set_visible(True) # riquadro nero
plt.show()
premere: "0" per inserire manualmente i parametri della malattia "1" per usare dei paramtri preimpostati 1
%matplotlib notebook
import networkx as nx
import EoN
import matplotlib.pyplot as plt
G = nx.grid_2d_graph(100,100)
initial_infections = [(u,v) for (u,v) in G if 45<u<55 and 45<v<55]
pos = {node:node for node in G}
sim_kwargs = {'pos': pos}
sim = EoN.fast_SIR(G, 2.0, 1.0, initial_infecteds = initial_infections,
tmax = 40, return_full_data=True, sim_kwargs = sim_kwargs)
ani=sim.animate(ts_plots=['I', 'SIR'], node_size = 4)
ani.save('SIR_2dgrid.mp4', fps=5, extra_args=['-vcodec', 'libx264'])
#fonte: https://epidemicsonnetworks.readthedocs.io/en/latest/Examples.html
import pygame, sys
import numpy as np
import random
BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE = (0, 100, 255)
RED = (255,0,0)
GREY = (230, 230, 230)
GREEN=(0,255,0)
BACKGROUND = BLACK
class Dot(pygame.sprite.Sprite):
def __init__(
self,
x,
y,
width,
height,
color=BLACK,
radius=5,
velocity=[0, 0],
randomize=False,
):
super().__init__()
self.image = pygame.Surface([radius * 2, radius * 2])
self.image.fill(BACKGROUND)
pygame.draw.circle(
self.image, color, (radius, radius), radius
)
self.rect = self.image.get_rect()
self.pos = np.array([x, y], dtype=np.float64)
self.vel = np.asarray(velocity, dtype=np.float64)
self.killswitch_on = False
self.recovered = False
self.randomize = randomize
self.WIDTH = width
self.HEIGHT = height
def update(self): # se un dot esce dal riquadro rientra dalla parte opposta
self.pos += self.vel
x, y = self.pos
if x < 0:
self.pos[0] = self.WIDTH
x = self.WIDTH
if x > self.WIDTH:
self.pos[0] = 0
x = 0
if y < 0:
self.pos[1] = self.HEIGHT
y = self.HEIGHT
if y > self.HEIGHT:
self.pos[1] = 0
y = 0
self.rect.x = x
self.rect.y = y
vel_norm = np.linalg.norm(self.vel)
if vel_norm > 3:
self.vel /= vel_norm
if self.randomize:
self.vel += np.random.rand(2) * 2 - 1 #-1/1 gfg
if self.killswitch_on:
self.cycles -= 1
if self.cycles<= 0:
self.killswitch_on = False # non si può morire 2 volte
some_number = np.random.rand()
if self.mortality_rate > some_number:
self.kill()
else:
self.recovered = True
def respawn(self, color, radius=5):
return Dot(
self.rect.x,
self.rect.y,
self.WIDTH,
self.HEIGHT,
color=color,
velocity=self.vel,
)
def killswitch(self, cycles=20, mortality_rate=0.2): # set di update di ciclo
self.killswitch_on = True
self.cycles = cycles
self.mortality_rate = mortality_rate
class Explosion(pygame.sprite.Sprite):
def __init__(self, x,y ):
pygame.sprite.Sprite.__init__(self)
self.x = x
self.y = y
self.images = []
for num in range(0,9):
img=pygame.image.load(f"regularExplosion0{num}.png")
img= pygame.transform.scale(img, (100,100))
self.images.append(img)
self.index = 0
self.image=self .images[self.index]
self.rect = self.image.get_rect()
self.rect.center = [x,y]
self.counter = 0
def update(self):
explosions_speed=4
self.counter +=1
if self.counter>=explosions_speed and self.index < len(self.images)-1:
self.counter=0
self.index+=1
self.image=self.images[self.index]
if self.index>= len(self.images) -1 and self.counter>=explosions_speed:
self.kill()
explosions_group = pygame.sprite.Group()
class Simulation:
def __init__(self, width=600, height=480):
self.WIDTH = width
self.HEIGHT = height
self.susceptible_container = pygame.sprite.Group()
self.infected_container = pygame.sprite.Group()
self.recovered_container = pygame.sprite.Group()
self.all_container = pygame.sprite.Group()
self.n_susceptible = 20
self.n_infected = 1
self.n_died=0
self.n_infected = 0
self.T = 1000
self.cycles = 20
self.mortality_rate = 0.2
def start(self, randomize=False):
self.N = (
self.n_susceptible + self.n_infected + self.n_quarantined
)
pygame.init()
screen = pygame.display.set_mode([self.WIDTH, self.HEIGHT])
for i in range(self.n_susceptible):
x = np.random.randint(0, self.WIDTH + 1)
y = np.random.randint(0, self.HEIGHT + 1)
vel = np.random.rand(2) * 2 - 1
guy = Dot(
x,
y,
self.WIDTH,
self.HEIGHT,
color=WHITE,
velocity=vel,
randomize=randomize,
)
self.susceptible_container.add(guy)
self.all_container.add(guy)
for i in range(self.n_quarantined):
x = np.random.randint(0, self.WIDTH + 1)
y = np.random.randint(0, self.HEIGHT + 1)
vel = [0, 0]
guy = Dot(
x,
y,
self.WIDTH,
self.HEIGHT,
color=BLUE,
velocity=vel,
randomize=False,
)
self.susceptible_container.add(guy)
self.all_container.add(guy)
for i in range(self.n_infected):
x = np.random.randint(0, self.WIDTH + 1)
y = np.random.randint(0, self.HEIGHT + 1)
vel = np.random.rand(2) * 2 - 1
guy = Dot(
x,
y,
self.WIDTH,
self.HEIGHT,
color=RED,
velocity=vel,
randomize=randomize,
)
self.infected_container.add(guy)
self.all_container.add(guy)
clock = pygame.time.Clock()
for i in range(self.T):
for event in pygame.event.get():
if event.type == pygame.QUIT:
sys.exit()
self.all_container.update()
screen.fill(BACKGROUND)
# Infetti
collision_group = pygame.sprite.groupcollide( # controllo collisione
self.susceptible_container,
self.infected_container,
True, #remove susc
False, #lascia inf
)
for guy in collision_group:
new_guy = guy.respawn(RED) #dobbiamo cambiare il colore -> rimuovo da S aggiungo a I
new_guy.vel *= -1
new_guy.killswitch(
self.cycles, self.mortality_rate
)
self.infected_container.add(new_guy)
self.all_container.add(new_guy)
# rimessi
recovered = []
for guy in self.infected_container:
some_number = np.random.rand()
if guy.recovered and self.mortality_rate > some_number:
new_guy = Explosion(random.randrange(0, 600),random.randrange(0, 480))
self.recovered_container.add(new_guy)
self.all_container.add(new_guy)
recovered.append(guy)
for guy in self.infected_container:
if guy.recovered:
new_guy = guy.respawn(GREEN)
self.recovered_container.add(new_guy)
self.all_container.add(new_guy)
recovered.append(guy)
if len(recovered) > 0:
self.infected_container.remove(*recovered) # leviamo i recovered da infetti
self.all_container.remove(*recovered)
self.all_container.draw(screen)
pygame.display.flip()
clock.tick(30)
pygame.quit()
if __name__ == "__main__":
scelto = int(input('premere:\n"0" per inserire manualmente i parametri della malattia\n"1" per usare dei paramteri preimpostati \nScelto: '))
if scelto == 0:
SIR = Simulation(600, 480)
SIR.n_susceptible = int(input('inserire il numero di persone sane: '))
SIR.n_quarantined = 0
SIR.n_infected = int(input('inserire il numero di persone infette: '))
SIR.cycles= 200 # indica sopo quanto tempo i dot possono iniziare a passare a R
SIR.mortality_rate = float(input('inserire il tasso di mortalità: '))
SIR.start(randomize=True)
elif scelto==1:
SIR = Simulation(600, 480)
SIR.n_susceptible = 100
SIR.n_quarantined = 0
SIR.n_infected = 5
SIR.cycles= 200
SIR.mortality_rate = 0.2
SIR.start(randomize=True)
else:
print('errore')
premere: "0" per inserire manualmente i parametri della malattia "1" per usare dei paramteri preimpostati Scelto: 1